{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DG求解一维Burgers方程"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DG单元操作"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dmatrix (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DifferentialEquations, Plots, LinearAlgebra, BenchmarkTools, Sundials\n",
    "\n",
    "function legendregausslobatto(m)\n",
    "    # Compute the m'th order Legendre-Gauss-Lobatto quadrature points, x, and weights, w\n",
    "    x = zeros(m+1); w = zeros(m+1)\n",
    "\n",
    "    if m==1\n",
    "        x[1]=-1.0; x[2]=1.0\n",
    "        w[1]= 1.0; w[2]=1.0\n",
    "        return x\n",
    "    end\n",
    "\n",
    "    if m==2\n",
    "        x[1]=-1.0; x[2]=0.0; x[3]=1.0\n",
    "        w[1]= 1/3; w[2]=4/3; w[3]=1/3\n",
    "        return x\n",
    "    end\n",
    "\n",
    "    J = zeros(m-1, m-1); h1 = 2*(0:m-2).+2\n",
    "    aux = broadcast(sqrt, (1:m-2).*((1:m-2).+ 2).*((1:m-2).+1).*((1:m-2).+1)./(h1[1:m-2].+1)./(h1[1:m-2].+3))\n",
    "    J[1:m-2,2:m-1] = Diagonal(2.0 ./(h1[1:m-2].+2).*aux)\n",
    "    J[1,1] = 0.0; J = J + J'\n",
    "\n",
    "    D, V = eigen(J)\n",
    "    i = sortperm(D); x = [-1.0; D[i]; 1.0]\n",
    "    x\n",
    "end\n",
    "\n",
    "function legendre(x,m)\n",
    "    # Evaluate orhonormal m'th order Legendre Polynomial at point x\n",
    "    xp = x\n",
    "    # Initial values P_0(x) and P_1(x)\n",
    "    PL = zeros(m+1, length(xp))\n",
    "    PL[1,:]  .= sqrt(1.0/2.0)\n",
    "    if  m == 0\n",
    "        P = PL'\n",
    "        return P\n",
    "    end\n",
    "    PL[2,:] = sqrt(3.0/2.0)*xp\n",
    "    if  m == 1\n",
    "        P = PL[m+1,:]\n",
    "        return P\n",
    "    end\n",
    "\n",
    "    # Forward recurrence using the symmetry of the recurrence.\n",
    "    aold = sqrt(1.0/3.0)\n",
    "    for i=1:m-1\n",
    "        anew = 2/(2*i+2)*sqrt((i+1)*(i+1)*(i+1)*(i+1)/(2*i+1)/(2*i+3))\n",
    "        PL[i+2,:] = 1/anew*(-aold*PL[i,:] + xp.*PL[i+1,:])\n",
    "        aold = anew\n",
    "    end\n",
    "    P = PL[m+1,:]\n",
    "    return P\n",
    "end\n",
    "\n",
    "function vandermonde(m,r)\n",
    "    # Initialize the 1D Vandermonde Matrix, V_{ij} = phi_j(r_i)\n",
    "    V = zeros(length(r),m+1)\n",
    "    for j in 1:m+1\n",
    "        V[:,j] = legendre(r, j-1)\n",
    "    end\n",
    "    V\n",
    "end\n",
    "\n",
    "function gradlegendre(r, m)\n",
    "    \"\"\"\n",
    "    Evaluate the derivative of the\n",
    "    m'th order Legendre polynomial at points r\n",
    "    \"\"\"\n",
    "    dP = zeros(length(r),1)\n",
    "    if m>0\n",
    "        Ph = -m*r.*legendre(r,m) + m*sqrt((2*m+1)/(2*m-1))*legendre(r,m-1)\n",
    "        dPe = r.^(m+1)*m*(m+1)/2*sqrt((2*m+1)/2)\n",
    "        endp = (broadcast(abs, broadcast(abs,r).-1)) .> 10*2.2204e-16\n",
    "        rh = r.*endp\n",
    "        dP = broadcast(~,endp).*dPe + endp.*Ph./(1.0 .- rh.^2)\n",
    "    end\n",
    "    dP\n",
    "end\n",
    "\n",
    "function dmatrix(m,r,V)\n",
    "    \"\"\"\n",
    "    Purpose : Initialize the (r) differentiation matrices,\n",
    "    evaluated at (r) at order m\n",
    "    \"\"\"\n",
    "    Vr = zeros(length(r),m+1)\n",
    "    for i=0:m\n",
    "       Vr[:,i+1] = gradlegendre(r,i)\n",
    "    end\n",
    "    D = Vr/V\n",
    "    D\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 右端项（Right Hand Side）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rhs (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rhs(u,p,t)\n",
    "    D, Φ, h, α, du, q, dq = p\n",
    "    \n",
    "    # Define u field differences at faces\n",
    "    du[1,2:end]     = (α/2.0 - 1)*(u[1,2:end]-u[end,1:end-1]) \n",
    "    du[end,1:end-1] = α/2.0*(u[end,1:end-1]-u[1,2:end])\n",
    "\n",
    "    # Impose boundary condition\n",
    "    uin  = - u[1,1]\n",
    "    uout = - u[end,end]\n",
    "    du[1,1] = (α/2.0 - 1)*(u[1,1] - uin)\n",
    "    du[end,end] = α/2.0*(u[end,end] - uout)\n",
    "    \n",
    "    # Compute q \n",
    "    q = -2.0/h*(D*u .- Φ*du) \n",
    "    \n",
    "    # Define q field differences at faces\n",
    "    dq[1,2:end]     = (α/2.0 - 1)*(q[1,2:end]-q[end,1:end-1]) \n",
    "    dq[end,1:end-1] = α/2.0*(q[end,1:end-1]-q[1,2:end])\n",
    "    \n",
    "    # Impose boundary condition\n",
    "    qin  = q[1,1]\n",
    "    qout = q[end,end]\n",
    "    dq[1,1] = (α/2.0 - 1)*(q[1,1] - qin)\n",
    "    dq[end,end] = α/2.0*(q[end,end] - qout)\n",
    "    \n",
    "    return -2.0/h*(D*q .- Φ*dq) \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题定义"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×100 Array{Float64,2}:\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0\n",
       " 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 2      # order of polynomials\n",
    "Np = p+1   # number of points\n",
    "K = 100    # number of elements\n",
    "xmin = 0\n",
    "xmax = 2π\n",
    "L = xmax - xmin\n",
    "h = L/K\n",
    "r = legendregausslobatto(p)\n",
    "VX = [xmin + i * h for i = 0:K]\n",
    "@inbounds x = ones(p + 1, 1) * VX[1:K]' + 0.5 * (r .+ 1) * (VX[2:K+1] - VX[1:K])'\n",
    "@inbounds V = vandermonde(p, r)\n",
    "@inbounds D = dmatrix(p, r, V)\n",
    "@inbounds Φ = V * V'\n",
    "@inbounds u0 = sin.(x) \n",
    "q  = zeros(Np, K)\n",
    "du = zeros(Np, K)\n",
    "dq = zeros(Np, K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 显式时间格式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  4.489733 seconds (8.43 M allocations: 1.660 GiB, 13.68% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.157025418563407e-5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tspan = (0.0, 0.8)\n",
    "prob = ODEProblem(rhs, u0, tspan, (D, Φ, h, 1, du, q, dq))\n",
    "@time uf = solve(prob,Tsit5(),save_everystep=false)\n",
    "error = norm(uf[2] - exp(-tspan[2])*sin.(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 隐式时间格式"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
